Stickleback example

Reading the shape information

Stickleback are swimbladdered fish. We have to shape files, one for the fish body and one for the swimbladder.
The shape file for the fish body contains one column named x_fb the fishbody (fb) along the x axis (Length), w_fb the width (seen from top) of the fish body at each position x (along the fish body length), z_fbL the lower height of the fish body along the x axis (lower (L) fish body (fb) extend in direction z) and z_fbU the upper height along the x axis (upper (U) fish body (fb) extend in direction z).

Analoguously, the swimbadder (sb) shape file contains x_sb - x axis (Length direction), w_sb - Width of the swimbladder along the x axis, z_sbL - lower swimbladder extend and z_sbU - upper swimbladder extend.

In our special case the csv files contain two different fish body and swimbladder shapes, identified by an additional ind column.

library(dplyr)

## 
## Attache Paket: 'dplyr'

## Die folgenden Objekte sind maskiert von 'package:stats':
## 
##     filter, lag

## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     intersect, setdiff, setequal, union

library(KRMr)

## Warning: vorhergehender Import 'ggplot2::last_plot' durch 'plotly::last_plot'
## während des Ladens von 'KRMr' ersetzt

## Warning: vorhergehender Import 'DT::dataTableOutput' durch
## 'shiny::dataTableOutput' während des Ladens von 'KRMr' ersetzt

## Warning: vorhergehender Import 'DT::renderDataTable' durch
## 'shiny::renderDataTable' während des Ladens von 'KRMr' ersetzt

stb_fb <- read.csv(paste0(dirname(getwd()),'/data/Sticklebacks_FB.csv'))
stb_sb <- read.csv(paste0(dirname(getwd()),'/data/Sticklebacks_SB.csv'))

library(plotly)

## Lade nötiges Paket: ggplot2

## 
## Attache Paket: 'plotly'

## Das folgende Objekt ist maskiert 'package:ggplot2':
## 
##     last_plot

## Das folgende Objekt ist maskiert 'package:stats':
## 
##     filter

## Das folgende Objekt ist maskiert 'package:graphics':
## 
##     layout

cyl <- function(r=c(1,2),r2=c(0.1,0.5),h=2,a=0, nt=10, nv=10){
  theta = seq(0,2*pi, length.out=nt)
  v = seq(a, a+h, length.out=nv)
  th = matrix(rep(theta,each=nv), ncol=nt)
  x = cos(th)*seq(r[1],r[2] , length.out=nv)
  y = sin(th)*seq(r2[1], r2[2] , length.out=nv)
  z = matrix(rep(v,times=nt), ncol=nt)
  return(list(x,y,z))
}

p = plot_ly()
fb0 = stb_fb[stb_fb$ind==58,]
sb0 = stb_sb[stb_sb$ind==58,]
for(i in seq(1, nrow(fb0)-1)){
  print(i)
  cc = cyl(r=c(fb0$w_fb[i], fb0$w_fb[i+1]), 
         a = fb0$x_fb[i],#abs(x_fb[i+1] - x_fb[i]), 
         h=fb0$x_fb[i+1]-fb0$x_fb[i],
         r2=c(fb0$z_fbU[i], fb0$z_fbU[i+1]))#x_fb[i])
  p = p%>%add_trace(x=cc[[3]], y=cc[[1]], z=cc[[2]], 
                    type='surface', 
                    showscale = FALSE, 
                    colorscale=list(c(0, 1), c("red", "red")),
                    opacity=0.5)
  
}

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8
## [1] 9
## [1] 10
## [1] 11
## [1] 12
## [1] 13
## [1] 14
## [1] 15
## [1] 16
## [1] 17
## [1] 18
## [1] 19
## [1] 20
## [1] 21
## [1] 22
## [1] 23
## [1] 24
## [1] 25
## [1] 26
## [1] 27
## [1] 28
## [1] 29
## [1] 30
## [1] 31
## [1] 32
## [1] 33
## [1] 34
## [1] 35
## [1] 36

for(i in seq(1, nrow(sb0)-1)){
  print(i)
  cc = cyl(r=c(sb0$w_sb[i], sb0$w_sb[i+1]), 
         a = sb0$x_sb[i],#abs(x_fb[i+1] - x_fb[i]), 
         h=sb0$x_sb[i+1]-sb0$x_sb[i],
         r2=c(sb0$z_sbU[i], sb0$z_sbU[i+1]))#x_fb[i])
  p = p%>%add_trace(x=cc[[3]], y=cc[[1]], z=cc[[2]], 
                    type='surface', 
                    showscale = FALSE, 
                    colorscale=list(c(0, 1), c("blue", "blue")),
                    opacity=0.5)
  
}

## [1] 1
## [1] 2
## [1] 3
## [1] 4

p%>%layout(scene=list(aspectmode='data'))

Let’s have a look at the shapes we loaded.
To plot the loaded shapes, we can use the shplot(x_fb, w_fb, x_sb, w_sb, z_fbU, z_fbL, z_sbU, z_sbL) function:

for (i in unique(stb_fb$ind)){
  fb =stb_fb%>%filter(ind==i)
  sb =stb_sb%>%filter(ind==i)
  #layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))
  par(mfrow=c(1,2))
  KRMr::shplot(x_fb = fb$x_fb, w_fb = fb$w_fb, 
             x_sb = sb$x_sb, w_sb = sb$w_sb,
             z_fbU = fb$z_fbU, z_fbL = fb$z_fbL,
             z_sbU = sb$z_sbU, z_sbL = sb$z_sbL)
}

For now, lets just include one shape example.

id = unique(stb_fb$ind)[1]
fb =stb_fb%>%filter(ind==i)
sb =stb_sb%>%filter(ind==i)

Running a KRM simulation

Now that we have our shapes sorted, let’s have a go at a KRM simulation.
We want to get TS at 38, 70, 120 and 200 kHz, with the below defined settings:

c.w = ambient water sound speed
rho.w = density ambient water
theta = orientation in degrees
c.fb = sound speed inside fish body
c.sb = sound speed inside swimbladder
rho.sb = density inside swimbladder
rho.fb = density inside fish body = L = length of the fishbody in m
x_fb = Fish body coordinates along x axis
x_sb = Swimbladder coordinates along x axis
w_fb = Fish body width along x axis
w_sb = Swimbladder width coordinates along x axis
z_fbU = Upper height of the fish body along the x axis
z_fbL = Lower height of the fish body along the x axis
z_fbU = Upper height of the swimbladder along the x axis
z_fbL = Lower height of the swimbladder along the x axis

A description of the parameters is also available in the help files ?KRMr::krm (runs a single krm simulation) or ?KRMr::krm.sim (runs multiple simulations if any input parameter contains more than one value (except for coordinates)

TS = KRMr::krm.sim(frequency =c(38,70,120,200) * 1000,
                   c.w = 1490,
                   rho.w = 1030,
                   theta=90,
                   c.fb = 1570,
                   c.sb = 345,
                   rho.sb = 1.24,
                   rho.fb = 1070,
                   L=0.25,
                   x_fb = fb$x_fb,
                   x_sb = sb$x_sb,
                   w_fb = fb$w_fb,
                   w_sb = sb$w_sb,
                   z_fbU = fb$z_fbU,
                   z_fbL = fb$z_fbL,
                   z_sbU = sb$z_sbU,
                   z_sbL = sb$z_sbL)

Let’s have a look at the outcome:

library(ggplot2)
ggplot(data=TS, 
       aes(x=frequency/1000, y=TS))+
  geom_line(size=1.2)+
  ylab(expression(TS~(~dB~re~m^2)))+
  xlab('Frequency (kHz)')+
  theme_classic()+
  theme(text=element_text(size=14),
        legend.position='top')

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

We can also compute slightly more sophisticated simulations of a range of frequencies, lengths and orientations (remermbe the KRM model is only valid for orientations, close to 90 degrees ~65-115 degrees):

TS = krm.sim(frequency =c(38,70,120,200) * 1000,
                   c.w = 1490,
                   rho.w = 1030,
                   theta=seq(65,115),
                   c.fb = 1570,
                   c.sb = 345,
                   rho.sb = 1.24,
                   rho.fb = 1070,
                   L=seq(0.01,0.1,by=0.01),
                   x_fb = fb$x_fb,
                   x_sb = sb$x_sb,
                   w_fb = fb$w_fb,
                   w_sb =  sb$w_sb,
                   z_fbU = fb$z_fbU,
                   z_fbL = fb$z_fbL,
                   z_sbU = sb$z_sbU,
                   z_sbL = sb$z_sbL)

Let’s plot the results:

ggplot2::ggplot(data=TS, ggplot2::aes(x= theta, y=L*1000 , z=TS))+
  facet_wrap(.~frequency/1000)+
  geom_contour_filled(bins=15)+
  ggplot2::scale_fill_viridis_d(expression(TS~'('~dB~re~ m^2~')'))+
  ggplot2::xlab('Theta (degrees)')+
  ggplot2::ylab(expression(L~'(mm)'))+
  ggplot2::scale_x_continuous(expand=c(0,0))+
  ggplot2::scale_y_continuous(expand=c(0,0))+
  ggplot2::theme_classic()+
  ggplot2::theme(legend.position='top',
                 legend.text=element_text(size=10),
                 text=element_text(size=16),
                 strip.background = element_blank(),
                 strip.text = element_text(face='bold'))

Let’s do the same without a swimbladder for a hypothetical fish without swimbladder (we simply assign the swimbladder coordinates NULL) for fish sizes of 1 cm to 10 cm without swimbladder:

TS = krm.sim(frequency =c(38,70,120,200) * 1000,
                   c.w = 1490,
                   rho.w = 1030,
                   theta=seq(65,115),
                   c.fb = 1570,
                   c.sb = 345,
                   rho.sb = 1.24,
                   rho.fb = 1070,
                   L=seq(0.01,0.1,by=0.01),
                   x_fb = fb$x_fb,
                   x_sb = NULL,
                   w_fb = fb$w_fb,
                   w_sb =  NULL,
                   z_fbU = fb$z_fbU,
                   z_fbL = fb$z_fbL,
                   z_sbU = NULL,
                   z_sbL = NULL)

Let’s plot the results:

ggplot2::ggplot(data=TS, ggplot2::aes(x= theta, y=L*1000 , z=TS))+
  facet_wrap(.~frequency/1000)+
  geom_contour_filled(bins=15)+
  ggplot2::scale_fill_viridis_d(expression(TS~'('~dB~re~ m^2~')'))+
  ggplot2::xlab('Theta (degrees)')+
  ggplot2::ylab(expression(L~'(mm)'))+
  ggplot2::scale_x_continuous(expand=c(0,0))+
  ggplot2::scale_y_continuous(expand=c(0,0))+
  ggplot2::theme_classic()+
  ggplot2::theme(legend.position='top',
                 legend.text=element_text(size=10),
                 text=element_text(size=16),
                 strip.background = element_blank(),
                 strip.text = element_text(face='bold'))

We can also look at orientation plots for a fish of 5 cm without a swimbladder at 38 kHz:

TSrot = krm.sim(frequency =c(38,200) * 1000,
                   c.w = 1490,
                   rho.w = 1030,
                   theta=65:115,
                   c.fb = 1570,
                   c.sb = 345,
                   rho.sb = 1.24,
                   rho.fb = 1070,
                   L=0.08,
                   x_fb = fb$x_fb,
                   x_sb = NA,#sb$x_sb,
                   w_fb = fb$w_fb,
                   w_sb = NA,#sb$w_sb,
                   z_fbU = fb$z_fbU,
                   z_fbL = fb$z_fbL,
                   z_sbU = NA,#sb$z_sbU,
                   z_sbL = NA)#sb$z_sbL)
ggplot(TSrot, aes(x = theta, y = TS, group=frequency/1000, col=TS)) +
  geom_path(size=1.2) +
  facet_wrap(.~frequency/1000)+
  scale_x_continuous(limits=c(0,360), breaks=c(0,90,180))+
  coord_polar(start=-pi/2,direction=1)+
  scale_colour_viridis_c(name='', limits=c(-80,-70), oob=scales::squish)+

  ylab(expression(TS~'('~dB~re~m^-2~')' ))+
  xlab(expression(paste(theta,' (\u00B0) ')))+
  
  geom_vline(xintercept=90, lty=2)+
  geom_vline(xintercept=0,size=1)+
  geom_vline(xintercept=180,size=1)+
  
  theme_classic()+
  theme(strip.background = element_blank(),
        strip.text=element_text(size=18),
        legend.position='top',
        legend.text=element_text(angle=-15),
        legend.key.width = unit(2,'cm'),
        axis.line.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.x = element_blank(),
        panel.spacing.y = unit(-5,'lines'),
        axis.title.x = element_text(vjust=27, size=16),
        axis.title.y = element_text(hjust=0.75, size=16),
        text=element_text(size=18))

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.